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We present a generalization of the coupled dipole method to the scattering of light by arbitrary 
periodic structures. This new formulation of the coupled dipole method relies on the same direct- 
space discretization scheme that is widely used to study the scattering of light by finite objects. 
Therefore, all the knowledge acquired previously for finite systems can be transposed to the study 
of periodic structures. 



I. INTRODUCTION 

In its original form, the coupled dipole method (CDM) 
was developed for the study, in free-space, of the scatter- 
ing of light by an object with finite dimensions. 1 ' 2 The 
method was subsequently extended to deal with objects 
near a substrate 3 ' 4 or inside a multilayer system. 5 The 
principle of the method is always the same: the object is 
represented by a cubic array of N polarizable subunits, 
each with a size small enough compared to the spatial 
variations of the electromagnetic field for the dipole ap- 
proximation to apply. If the CDM could be extended 
to deal with local scatterers near periodic structures, the 
CDM could then also be used, for example, to study light 
scattering by objects near surface gratings or by defects 
or cavities in photonic crystals. The first step toward 
such an extension is to develop a form of the CDM ca- 
pable of describing periodic structures efficiently. In this 
paper, we present a generalization of the CDM to arbi- 
trary periodic structures. 



II. SELF-CONSISTENT FIELD FOR A PERIODIC 
STRUCTURE 

We consider a plane substrate occupying the region 
z < 0. For a single object on the substrate, the self- 
consistent field at the i th subunit at location is given 
by 

N 

E(r;,w) = E (r. i ,w) +^[S(r. i ,r J ,w) 
j'=i 

+ F(r i ,r j ,u))]a j (u>)E(r j ,u). (1) 

where Eo(rj,cj) is the (initial) field at in the absence 
of the scattering object. Note that none of the subunits 
lies in the plane z = 0. The tensors F and S are the 



field susceptibilities (linear responses) associated with 
the free space 6 and the substrate. 7 on(pS) is the dynamic 
polarizability of the i th subunit and includes radiation 
reaction. 2 ' 8 The self-consistent field E(rj,a;) is found by 
solving the symmetric linear system formed by writing 
Eq. (1) for i = l,N. The total field at position r is 
computed as 



E(r, u) = E (r, uj) 



N 

]T[S(r,r,, W ) 



+ F(r,r i ,w)]a i (u;)E(r i ,w). (2) 

This conventional form of the CDM is well adapted to 
deal with localized objects. If, instead of a single object, 
one wants to study a periodic structure created by the 
repetition of the object over a lattice located above the 
substrate, Eq. (1) becomes 

E(r,,cj) = E (Ti,u) 

N oo 

+ [S(rj,r,- +mu + nv,u>) 

j — l ra.n— — oc 



F(r, 



iu + nv, w)]cvj(w)E(rj + mu + nv, uj). (3) 



Eykp|| 




FIG. 1. Example of a periodic structure created by the 
repetition of an object over a lattice parallel to a substrate. 
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The vectors u and v are the basis vectors of the lattice 
(Fig. 1). The index i runs over all the subunits of the 
structure, is the position of subunit i. The sum over j 
is restricted to the N subunits of a single object with po- 
sition fj inside the object. The number of subunits is now 
infinite, and therefore so is the size of the linear system to 
be solved. One solution would be to truncate the infinite 
sum and solve the system for a large but finite number 
of objects, but this is impractical because the sums over 
the lattice converge very slowly. This problem can be cir- 
cumvented by using a plane-wave decomposition of the 
incident field. In the case of plane-wave (propagating or 
evanescent) illumination, the field above the surface can 
be written as (we note by k || the projection of vector k 
on a plane parallel to the surface) 



/oo 
dr \\ H 5(r|| -mu-nv)exp(iko||.r||) 

Til -Tl — — OO 

x [S(p ij -r^z i ,z j ,u))+F(p ij -r^,z i ,z j ,u/)] (6) 

We define the two-dimensional Fourier transform as : 
•^*[fr( r ||)] — / ^ r ||K r ||) ex P(~ * r ||-h||) 5 an( i i ts inverse as 
^-^(h,!)] = 1/(2tt) 2 x /dh||B(h||)exp(ir||.h||). Using 
the Parseval-Plancherel theorem Eq.(6) becomes 

If °° 
K = T^J dh U M E <5(h||-mu'-nv' + k ||) 

m.n— — 00 

x T[S{pij - ru,z i ,z j ,u>)+F(p ij - r«,Zi,z h u)\ (7) 



where u' = 2ir(v y yL — v x y)/(u x v y 



c u y ) and v' 



xti z_ \ m /- \ r-i / m 27r(— u v & + u x y)/(u r v v — v r u v ) are the basis vectors of 

E (r j +mu + nv,o;)=Eo(r i , W )exp[ t ko||.(mu + nv)] ^ J f"\* » L v lr^2 lu , _ „, .. , ~ 



(4) 

where ko is the wave vector in free space. Because of the 
periodicity of the system and the translational invariance 
of the field susceptibilities, the self consistent field satis- 
fies the same relation as the incident field (Eq. (4)), and 
at any subunit Eq. (3) can be written as 

E(r i5 w) = E (r 4 ,w) 

N / 00 

+ X! [S(Ti,fj+mu + nv,w) 

j—l \ m,n=— 00 



+ F(rj, Vj + rau + nv, ui)] cxp[ik || .(mu + nv)] 
ajME(Fj,w). 



(5) 



The self-consistent field on the right-hand side of Eq.(5) 
is independent of (m,n) and can be taken out of the in- 
finite sum. Hence the sum over subunits in Eq.(5) only 
involves j = 1, N, that is the number of subunits in a unit 
cell, which we choose to be the cell for which m = n = 0. 
Moreover, because of the translational symmetry of the 
self-consistent field, we only need to find E in one cell. 
Once the self-consistent field is found in the central cell, 
the field in any other cell is obtained by multiplying by 
the appropriate phase factor. Thus we only have to solve 
a linear system of the same size as the one describing a 
single object. The major issue in solving Eq. (5) is to 
compute efficiently the infinite, slowly convergent sums 
without performing a truncation of the sums. This is 
possible owing to the translational invariance of the field- 
susceptibilities in a plane parallel to the surface. The 
dependence on (Fj, Fj, ui) can be written as (pij, Zi, Zj, ui) 
with p^ = (Fj — Fj)||. Hence, the infinite sums of Eq.(5) 
become: 

00 

K = ^ [S(Fj,Fj +mu + nv, ui) 

m.n— — 00 

+ F(Fj, fj + mu + nv, u))] cxp[ik || .(mu + nv)] 



the reciprocal lattice, and M = (2n) / (u x v y — v x u y ). x 
and y are the basis vectors of the coordinate system. Us- 
ing the angular spectrum representations W and G of 
tensors S and F Eq.(7) becomes 6,7 



K = ^-M ^2 exp[z(mu' + nv' + k ||).py] x 



m.n— — 00 



{W(mu' + nv' + k || , k ) exp[iw {zi + Zj)] 
+ G(mu' + nv' + k ||,k ) exp[iw \zi - Zj\]}, 



with 



G(k|,,ko) 



Wo 



k 2 W k 2 



Wo W(, 



) 



-l k v 



(8) 



(9) 



and 



W(k||,k ) = 

' ' k 2 w 



k,.k 



wok' 



% (w^A p + k 2 A s ) k x A p \ 



ktkiA s 



k 2 .w 



-k A 



-k A 



kyAp 



(10) 



where 7 = sign(z; — zj), k|| = mu'+nv'+k || = k x St+k y y, 
and wo is the component along z of the wave vector ko , 
i.e., w = (k 2 - kl - ky) 1 / 2 . A p and A s are the Fresnel 
reflection coefficients for the substrate. Sums involving 
different susceptibility tensors (free-space or surface) will 
have a different behavior, due to the different arguments 
of the exponential terms (zi + Zj and | Zi — Zj \ ) . They will 
be computed separately. 

For the surface term, the convergence of the sum is en- 
sured by the exponential term. As m and n increase, the 
magnitude of kii increases and the nature of the plane 
wave changes from propagating to evanescent. Because 
Zi + Zj never vanishes, and because the subunits are never 
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exactly on the surface, this exponential term is always 
present and ensures the rapid convergence of the sum. 

For the free-space part, the argument of the exponen- 
tial term is \zi — Zj\ and the rapid convergence of the sums 
is not as trivial. We use the method introduced to derive 
the Green function of a 2D square grating. 9 We consider 
two cases. The first case pertains to the interaction be- 
tween elements from different "layers" of the lattice, and 
corresponds to the case Zi ^ Zj. This case is similar to 
the surface problem and the convergence of the sum is 
ensured by the exponential term. 

In the second and the exponential term 

disappears. We cast the free-space part of the infinite 
sum in two different forms. We note by a(ry , z% — Zj) the 
sum in the direct space (F terms in Eq.(6)). We note by 
A(k||, Zi — zj) the sum in the reciprocal space (G terms 
in Eq.(8)). When Zi — Zj we write the sum as 

a(r,| , 0) = A(k,j , h) + [a(r,| , 0) - a(r,| , h)], (11) 

where h is an offset parameter. The auxiliary sum 
A(k||,ft) can be computed efficiently owing to the pres- 
ence of an exponentially decreasing term. The difference 
of direct-space sums a(r||,0) — a(r^,h) goes as l/r| and 
can also be computed efficiently. With Eq.(ll) we can 
ensure a rapid convergence of the sums in a discretization 
plane despite the absence of an exponentially decreasing 
term in the original sum. 

To improve further the convergence of the sums, we use 
Shanks' accelerator. 10 Because we have two sums (over 
m and n) one solution would be to apply successively 
Shanks' accelerator to the inner (n) and outer (m) sums 
(as suggested in Ref. 11). The problem with this ap- 
proach is that in our case, the convergence of the inner 
sum (over n) can be very slow for high values of m (outer 
sum). A better solution consists in defining one element 
I of the Shanks's series as the sum over m — —I, I for 
n = —I, and n = —I, I for m = —1 + 1, I — 1. 
This strategy gets rid of the inner/outer sum problem 
and results in a faster convergence and an easier imple- 
mentation of the Shanks algorithm. 

Note that there is another way of computing efficiently 
the free-space term. As we did earlier, when we in- 
troduced a parameter h, it is possible to split the in- 
finite sum (F) terms in Eq.(6) in two parts; one in 
the direct space and one in the reciprocal space, where 
these two sums converge quickly owing to a damping 
fu nction. 12 ' 13 The convergence is the best when h = 
\Jit/{u x v y — v x u y ). Poppe et al. introduced this method 
to study the optical response of an atomic monolayer; the 
period of the structure was therefore very small compared 
to the wavelength. 

Once the periodic susceptibility tensors are known, 
we solve the linear system of Eq. (5) to find the self- 
consistent field at each site. Once the field at all subunits 
is known, the scattered field at any position r, above, be- 
low or inside the periodic structures is readily computed 
through Eq.(5) with the exchange r <-» r*. Notice that 



the new linear system is no longer symmetric. This is 
due to the fact that the elements of the system depend 
on the incident plane wave via the exponential term in 
Eq.(5). 



III. EXAMPLE: SCATTERING BY A PERIODIC 
STRUCTURE LYING ON A SUBSTRATE 

To illustrate the method we consider the case of a di- 
electric substrate (the relative permittivity is 2.25) on 
which lies a 2D grating of parallelepipeds with the same 
permittivity. The structure is illuminated in TM po- 
larization from the substrate side by total internal re- 
flection at an angle of incidence 8 — 45°; then k ii = 

(^sin6»V2^25,0). The wavelength in vacuum is A = 
632.8 nm, and the basis vectors of the lattice u = (a, 0), 
v = (0,a). 



5.2 




x (nm) 

FIG. 2. Intensity of the electric field above a dielectric sub- 
strate in the direction of the x-axis with a 2D grating of par- 
allelepipeds. The inset shows the geometry used. The solid 
line is for an isolated parallelepiped. The other curves are 
obtained for the 2D grating with a = 100 nm (dotted line), 
a — 200 nm (dot dashed line), a = 1000 nm (dashed line). 

The parallelepipeds have a square base of 40 x 40 nm 2 , 
and a height of 20 nm (see inset Fig. 2). In Fig. 2 we 
present the intensity of the electric field along the £-axis, 
60 nm above the dielectric substrate for different value of 
a. The curves are obtained for N = 256, hence the size of 
the subunit is 5 x 5 x 5 nm 3 (but convergence is already 
achieved for N = 32). Notice that the solid line is for an 
isolated parallelepiped on the substrate, i.e., the electric 
field is computed with the conventional CDM. 4 When a 
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is small, the computed curves for the electric field are no- 
tably different from the single object case. This denotes 
a strong coupling between parallelepipeds. Conversely 
for large a (a = 1000 nm), the curve is very similar to 
the curve for an isolated parallelepiped. 

Table 1 presents the computation time for the coeffi- 



cients of the linear system (Eq. (6)) used to solve Eq. 
(5), for different values of N, and three values of a. The 
factor h has an important influence on the computation 
time, therefore we have chosen the optimal value of h for 
each case. As a reference we use the conventional CDM 
to compute the field for a single parallelepiped. 14 



N 


32 


256 


500 


1372 


CDM 


2 


18 


39 


137 


a = 
100 nm 


CDMi 


0.3 (2) 


4(17) 


10 (34) 


43 (116) 


CDM 2 


0.2 (0.4) 


2.7 (5.5) 


7(13) 


29 (54) 


a = 
200 nm 


CDMi 


0.7 (2) 


12 (30) 


30 (75) 


119 (300) 


CDM 2 


0.4 (1) 


7(16) 


18 (40) 


72 (158) 


a = 
1 fj,m 


CDMi 


5.7 (16) 


96 (281) 


246 (684) 


949 (4020) 


CDM 2 


5.6 (16) 


96 (276) 


233 (674) 


900 (2460) 



TABLE I. Computation time in seconds for the coefficients of the linear system (Eq. (6)) used to solve Eq. (5). N is the 
number of subunits. CDM is the time for the classical CDM for one parallelepiped. CDMi is the time for the periodic CDM 
when the free space contribution is computed with Eq.(8), and CDM 2 is the time for the periodic CDM when the free space 
contribution is computed with Ref.[12]. The infinite sums of the series are stopped when the relative error is less than 10~ 3 

(io- 6 ). 



Table 1 shows three computation times: CDM is the 
time for the classical CDM for one parallelepiped. CDMi 
is the time for the periodic CDM when the free space 
contribution is computed with Eq.(8), and CDM 2 is the 
time for the periodic CDM when the free space contri- 
bution is computed with Rcf.[12]. CDM2 is faster than 
CDMi for small periods. For a — 1/zm the computa- 
tion times are similar. For larger periods CDM 2 fails to 
converge to the reference result because the method of 
Rcf. [12] used to compute the free space term does not 
work well for large a. We note that the computation 
time increases with a. This is due mainly to the surface 
term. The convergence of the series depends on the term 
exp[wo(zi + Zj)]. In our case the modulus of the vectors 
of the reciprocal basis are |u'| = |v'| = 27r/a. Hence when 
a decreases, the modulus of the vector basis increases, wq 
becomes imaginary for smaller values of (m, n), and the 
exponential term produces a stronger damping. Obvi- 
ously, when N increases, the computation time increases 
due to the increased number of subunits involved. But 
there is another effect of the surface term. As the size 
of the subunit becomes smaller (N increases), there are 
more subunits close to the substrate with a small value 
of Zi + Zj and a slower exponential decay. When we com- 
pare the classical CDM to the periodic CDM we see that 
for a smaller than 200 nm the computation time of the 
periodic CDM is shorter. When the size of the period be- 
comes larger than the wavelength used, the convergence 
becomes slower. 



IV. CONCLUSION 



In conclusion we have generalized the coupled dipole 
method (CDM) to periodic structures. We have dis- 
cussed explicitly the case of a three-dimensional struc- 
ture, periodic in two directions, placed on a substrate. 
However, the principle of the approach described here ap- 
plies to a broad range of configurations with one, two or 
three-dimensional structures. The main advantage of this 
new formulation is that it relics on the same straightfor- 
ward, direct-space discretization scheme that is used for 
a single localized object. Therefore, all the knowledge ac- 
quired previously in CDM modeling of finite systems can 
be transposed to the study of periodic structures. 15 Opti- 
cal anisotropy, for instance, can be included by taking the 
appropriate permittivity tensor. Also, as shown here, the 
symmetry of the periodic lattice can be arbitrary. Here, 
we have considered the case of plane wave illumination. 
In the case of arbitrary illumination, each spectral com- 
ponent of the incident field must be treated individually. 
An interesting extension of the present work would be 
to merge the periodic CDM and the conventional CDM 
into a single approach to light scattering. This would be 
particularly useful in dealing with localized defects in pe- 
riodic structures or the interaction between a near-field 
probe (microscope tip, fluorescing particle,...) and a pe- 
riodic system. The periodic generalization of the coupled 
dipole method can also be used to draw a better physical 
picture of the local-field corrections that appear during 
the multiple scattering of light by a discrete set of scat- 
terers. 16 

P. C. Chaumct's email address is pchaumet@loe.u- 
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